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' We describe new planetesimal accretion calculations in the Kuiper Belt that 

include fragmentation and velocity evolution. All models produce two power law 
cumulative size distributions, Nc oc r~ 2 " 5 for radii ^ 0.3-3 km and Nc oc r~ 3 for 
radii ^ 1-3 km. The power law indices are nearly independent of the initial mass 
in the annulus, M ; the initial eccentricity of the planetesimal swarm, eo; and 
the initial size distribution of the planetesimal swarm. The transition between 
the two power laws moves to larger radii as eo increases. The maximum size of 
objects depends on their intrinsic tensile strength, So; Pluto formation requires 
S Q ^ 300 erg g -1 . The timescale to produce Pluto-sized objects, t p , is roughly 
proportional to M _1 and eo, and is less sensitive to other input parameters. Our 
models yield Tp w 30-40 Myr for planetesimals with eo = 10 -3 in a Minimum 
Mass Solar Nebula. The production of several 'Plutos' and ~ 10 5 50 km radius 
Kuiper Belt objects leaves most of the initial mass in 0.1-10 km radius objects 
that can be collisionally depleted over the age of the solar system. These results 
resolve the puzzle of large Kuiper Belt objects in a small mass Kuiper Belt. 
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1. INTRODUCTION 



Recent discoveries of slow- moving objects beyond the orbit of Neptune have radically 
changed our understanding of the outer solar system. These observations have revealed a 
large population of Kuiper Belt objects (KBOs) in orbits with semi-major axes of 39-45 
AU QJewitt fc Luu 19931 , 1995; [Irwin et al. 1~995[ ; [Williams et al. 1995| ; Hewitt et al. 1996| ; 
Luu et al. 1997j |Gladman fc Kavelaars 1997] ; |Gladman et al. 1998|) . KBOs with reliable 



orbits have a cumulative size distribution that follows Nc oc r 
al. 1998| ; |Luu fc Jewitt 199~8| ). The estimated population of ^ 
km indicates a total mass, Mr bo ~ O.IM^Q, for reasonable assumptions about the albedo 
and distance distribution. 

The small mass in KBOs is a problem for current planet formation theories. In most 
theories, small planetesimals in the solar nebula grow by collisional accumulation (e.g., 
Safronov 1969| ; |Greenberg et al. 1978| , 1984; |Wetherill fc Stewart 1989] 1993; [Spaute et a* 



q , with q = 3 ± 0.5 QJewitt ei 
10 5 KBOs with radii r ^ 50 



1991|) . This growth eventually produces one or more 'cores' that accumulate most, if not 



all, of the solid mass in an annular 'feeding zone' defined by balancing the gravity of the 
growing planetesimal with the gravity of the Sun and the rest of the disk. Large cores with 
masses of 1-10 Me can also accrete gas from the feeding zone (Pollack 1984j) . Applied 
to the inner solar system, this model generally accounts for the masses of the terrestrial 
and gas giant planets (e.g., [Pollack et al. 199"6| , |Weidenschilling et al. 1997| but see |Bossj 
1997|) . At 30-50 AU, however, the timescale to produce planet-sized, ~ 1000 km, objects 
exceeds the disk lifetime, ^ 100 Myr, unless the mass of the outer disk is ~ 10 2 -10 3 Mkbo 



(e.g., [Fernandez fc lp 1981| , 1984; |lp 1989| ; |Stern 1995| , 1996; |Stern fc Colwell 1997a|) . The 
presence of large KBOs in a small mass Kuiper Belt is thus a mystery. 

In Kenyon & Luu (1998, hereafter KL98), we began to address this issue by considering 
planetesimal growth in a single annulus at 35 AU. We showed that calculations including 
accretion and velocity evolution naturally produce several "Plutos" with radii exceeding 
1000 km and numerous 50 km radius KBOs on timescales of r ~ 20-200 Myr for a wide 
range of initial conditions in plausible solar nebulae. These timescales indicate that Pluto 
can form in the outer solar system in parallel with the condensation of the outermost large 
planets. 

In this paper, we extend KL98 by adding fragmentation to our planetesimal evolution 
code. The code generally matches published calculations at 1 AU (Wetherill & Stewart 
1993, hereafter WS93) and at 40 AU (Davis & Farinella 1997). Our numerical results 
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demonstrate that fragmentation and velocity evolution damp runaway growth to provide a 
self-limiting mechanism for the formation of KBOs and Pluto in the Kuiper Belt. These 
calculations produce several Plutos and at least 10 5 50 km radius KBOs on timescales of 
2-100 Myr in annuli with modest surface densities, 2.0-0.14 g cm -2 , of solid material at 
35 AU. Our analysis sets a lower limit on the intrinsic strength of KBOs, ~ 300 erg g -1 , 
and indicates that the initial size distribution, the initial eccentricity of the planetesimal 
swarm, and the details of the fragmentation algorithm have little impact on the resulting 
size distribution of KBOs. 

Our results appear to resolve the mystery of large KBOs in a small mass Kuiper Belt. 
Planetesimal evolution at 35-50 AU in a Minimum Mass Solar Nebula ( jHayashi 198 1| ) 



naturally produces large KBOs in numbers close to those currently observed. Most of 
the disk mass ends up in smaller objects, 0.1-10 km, that can be collisionally depleted 
over the age of the solar system. This depletion rate depends on the intrinsic strength 
and eccentricity distribution of KBOs ([Stern 1996| ; |Davis fc Farinella 1997| ; [Stern fc| 



|Colwell 19974 1997b). Future observations can place better constraints on these physical 
parameters and provide additional tests of our interpretation of KBO formation. 

We outline the fragmentation model and tests in Sec. 2, describe calculations for 
the Kuiper Belt in Sec. 3, and conclude with a discussion and summary in Sec. 4. The 
Appendix contains a complete description of the fragmentation algorithms and updates of 
the coagulation code from KL98. 



2. The Model 

As in KL98, we adopt Safronov's (1969) particle- in- a-box method, where planetesimals 
are treated as a statistical ensemble of masses with a distribution of horizontal and 
vertical velocities about a Keplerian orbit. Our calculations begin with a differential mass 
distribution n{nii) in a single accumulation zone centered at a heliocentric distance a with 
an inner radius at a — Aa/2 and an outer radius at a + Aa/2. The horizontal and vertical 
velocity dispersions of these bodies are hi(t) and Vi(t). We approximate the continuous 
distribution of particle masses with discrete batches having particle populations rij(t) and 
total masses Mj(i) (WS93). The average mass of a batch, rriiit) = M^fj/n^t), changes 
with time as collisions add and remove bodies from the batch. This procedure naturally 
conserves mass and allows a coarser grid than calculations with fixed mass bins ( |WetheriT| 



19901 IWS93D . 



To evolve the mass and velocity distributions in time, we solve the coagulation and 



3 



energy conservation equations for an ensemble of objects with masses ranging from ~ 10 7 g 
to ~ 10 26 g. The Appendix of KL98 describes our treatment of accretion and velocity 
evolution, and compares numerical results with analytic solutions for standard test cases. In 
this paper, we add fragmentation to the problem. We adopt a simple treatment of collision 
outcomes, following Greenberg et al. (1978; see also Davis et al. 1985| , 1994; WS93): 



1. Mergers - two bodies collide and merge into a single object with no debris; 

2. Cratering - two bodies merge into a single object but produce debris with a mass much 
smaller than the mass of the merged object; 

3. Rebounds - two bodies collide and produce some debris but do not merge into a single 
object; and 

4. Disruption - two bodies collide and produce debris with a mass comparable to the mass 
of the two initial bodies. 



We consider two algorithms for treating the cratering and disruption of planetesimals. 
Both the WS93 and the Davis et al. (1984) algorithms estimate the amount of debris 
rrifjj produced from impacts with velocities exceeding a threshold velocity Vf. In general, 
rrifjj scales with the impact energy. WS93 assume that the debris has the same relative 
velocity V%j of the colliding bodies; Davis et al. (1984) assume that the debris receives a 
fixed fraction f KE of the impact kinetic energy. In both cases, we adopt a coefficient of 
restitution cr to allow for rebound collisions that produce debris but no single merged 
body. We follow Greenberg et al. (1978) and adopt separate values of cr for collisions with 
(c-z) and without (ci) cratering. The Appendix describes these algorithms in more detail. 

To test our fragmentation procedures, we attempt to duplicate WS93's calculations of 
planetary embryo formation at 1 AU. Their model begins with 8.33 x 10 8 planetesimals 
having radii of 8 km and a velocity dispersion of 4.7 m s _1 relative to a Keplerian orbit 
(Table 1; see also Table 1 of |WS93| ). Table 2 summarizes our results using the WS93 initial 
conditions with mass spacing factors of 5 = mj+i/mj = 1.12, 1.25, 1.4, and 2.0 between 
successive mass batches. We adopt the analytic cross-sections of Spaute et al. (1991), which 
yield results identical to the numerical cross-sections of WS93 (|KL98| ). Figure 1 shows our 



reproduction of WS93 for S = 1.25. This model produces twelve 3-7 x 10 25 g objects with 
velocity dispersions of 60-100 m s -1 in 1.15 x 10 5 yr. In contrast, WS93 produce seven 
1-3 x 10 26 g objects (see Figures 1-3 in WS93). Despite this factor of two difference in 
the total mass contained in the most massive bodies, our calculation has the same broad 
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"plateau" in the cumulative number distribution Nq at log mj = 24-26 and a similar power 
law dependence, N c oc mj\ at log rrii = 21-23. Our "fragmentation tail" at log = 7-18 
has the standard power law dependence, N c oc m~ 0,8 ( |Dohnanyi 1969| ), at large masses but 
flattens out more than the WS93 result at small masses. 

The evolution of particle velocities in our calculations generally agrees with WS93 
(Figure 2; see Figures 2-3 of WS93). All velocities increase monotonically with time due to 
viscous stirring. The velocities of large bodies grow very slowly, because dynamical friction 
transfers their kinetic energy to small bodies. The model maintains a nearly constant ratio 
of vertical to horizontal velocity, Vi/hi ~ 0.53, for all but the most massive bodies, which 
have Vi/hi < 0.5 (Figure 2, right panel). This result agrees with previous calculations 
(Barge & Pellat 1990, 1991; [Hornung et ai. 1985| ). At late times, our velocities for small 



bodies, hi w 200 m s 1 at m, ~ 10 7 g, are roughly a factor of two larger than those of 
WS93. Velocities for large bodies, hi ~ 100 m s" 1 at m ; ~ 10 25 g, roughly equal those of 
WS93. 

The higher velocities of our calculation lead to additional mass loss from gas drag and 
fragmentation. Gas drag removes material from the annulus; fragmentation produces bodies 
with masses less than the minimum mass, ~ 10 7 g, of our numerical grid. We typically 
lose ~ 25% of the initial mass to cratering and catastrophic fragmentation and another 
20%-25% to gas drag. WS93 lost a comparable amount of mass to fragmentation but only 
~ 5% to gas drag. 

The agreement between our results and those of WS93 depends on the mass spacing 
factor 5. For r <^ 30,000 yr, large 5 models take longer to produce 1000-2000 km objects 
than small 5 models. The 'lag' relative to our 5 = 1.12 model increases from l%-2% 
for 5 = 1.25 up to ~ 10% for 5 = 2 (see also |Uhtsuki fc Makagawa 198"E|; phtsuki et at. 



nM?t |Wetherill 199H| ; |Kolvoord k, Greenberg Iggg] KL98). The poor resolution of 8 > 1.4 



models delays the production of massive bodies that 'runaway' from the rest of the mass 
distribution. In WS93, most runaway bodies are also 'isolated' bodies that do not interact 
with one another but do interact with the rest of the mass distribution (see WS93 and 
the Appendix). Delays in the production of isolated, runaway bodies lead to 10%-20% 
increases in the velocity dispersion of all bodies. Larger velocities reduce gravitational 
focusing factors and slow the growth of the largest bodies. The mass of the largest object 
thus increases from 7-8 x 10 25 g for 5 = 1.25-1.4 to 9.5 x 10 25 g for 5 = 1.12. This trend 
suggests that calculations with 5 = 1.01-1.10, as in WS93, would improve the agreement 
between our results and those of WS93. 

The coagulation results also depend on the treatment of low velocity collisions. 
Our standard expressions for long-range gravitational interactions fail when the velocity 
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dispersion is close to the Hill velocity, vh ( [Ida 19901 ; [Barge fc Pellat 1991| ). We use Ida's 



(1990) n-body calculations and adopt his simple expressions to derive the velocity evolution 
for collisions in the low velocity limit, Vy < Vi v vh, where Vi v = 2-5 is a constant (see the 
Appendix). This change slows down the velocity evolution, because it introduces a lower 
limit to the timescales for dynamical friction and the inclination component of viscous 
stirring (|Ida 1990Q . We adopt V lv = 2, following WS93, for 1 AU calculations. Unlike WS93, 



however, the small bodies in our calculation do not quite reach the low velocity limit and 
increase in velocity throughout the evolution. A modest increase in our low velocity limit, 
V\ v = 3.5, leads to better agreement between the velocity behavior of our models and that 
of WS93. 

To test our procedures further, we repeat the 1 AU calculations using the Davis et al. 
(1984) fragmentation algorithm. Table 3 lists our results for J'ke = 0.1 and 5 = 1.25, 1.4, 
and 2.0. These calculations yield solutions similar to our WS93 models. Calculations with 
smaller 5 produce larger bodies at earlier times than models with large S (Table 3). The 
mass of the largest object increases from rrii ~ 7 x 10 25 g for 5 = 2.0 to wij « 1.2 x 10 26 
g for 5 = 1.4 to m, w 1.8 x 10 26 g for 5 = 1.25. The velocity dispersions of the small 
bodies are a factor of 2-3 larger at r ^ 30,000 yr because they receive a larger fraction of 
the fragmentation energy than in the WS93 algorithm. This trend reverses for r ^ 50,000 
yr, because the Davis et al. algorithm produces more debris for collisions between large 
objects than the WS93 algorithm. These large objects have small velocity dispersions, so 
the velocity dispersions of the small bodies increase more slowly than in the WS93 model. 
The mass lost from fragmentation, ~ 5% of the initial mass, and gas drag, ~ 20% of the 
initial mass, is also smaller. 

We conclude that our coagulation code with fragmentation and velocity evolution 
reproduces 'standard' calculations. The results in Tables 2-3 generally agree with published 
results, despite some differences in the evolution of the velocity dispersion and the radius 
of the largest object. These discrepancies probably result from subtle differences in the 
implementation of the algorithms. None seem significant, because the major results do not 
depend on the input parameters (see also WS93): all calculations produce several objects 
with rrii ~ 10 26 g that contain most of the remaining mass in the annulus. 
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3. Kuiper Belt Calculations 



3.1. Starting Conditions 

As in KL98, we rely on observations of other stellar systems and solar nebula 
models to choose appropriate initial conditions for our Kuiper Belt calculations. Recent 
observations indicate lifetimes of ~ 5-10 Myr for typical gaseous disks surrounding nearby 
pre-main-sequence stars and for the solar nebula QSargent fc Beckwith 1993"t |Russell e\ 
al. 1996| ; |Hartmann et al. 1998| ). In our solar system, Neptune formation places another 



constraint on the KBO growth time, because Neptune excites KBOs through gravitational 
perturbations. Recent calculations suggest Neptune can form in 5-100 Myr (|lp 1989 ; 



[Lissauer et al. 1996] ; |Pollack et al. 199"E| ). Once formed, Neptune inhibits KBO formation at 



30-40 AU by increasing particle random velocities on timescales of 20-100 Myr ( |Holman~S| 
Wisdom 1993] puncan et al. 1995| ; |Morbidelli fc Valsecchi I9"9T| ). We thus adopt 100 Myr 



as an upper limit to the KBO formation timescale at 30-40 AU. 

Our starting conditions for KBO calculations are similar to KL98. We consider an 
annulus centered at 35 AU with a width of 6 AU. This annulus can accommodate at least 
10-100 isolated bodies with rrii ^ 10 24 g for e < 0.01. Instead of the single starting radius 
used in KL98, the present calculations begin with No bodies in a size distribution with 
radii, = 1-80 m. This initial population has a power law cumulative size distribution, 
N c oc r~ go . We usually adopt q = 3; the final size distribution appears to depend very little 
on q , as discussed below. The planetesimals have a small initial eccentricity (Malhotra 



1995) that is independent of size and an equilibrium ratio of inclination to eccentricity, 
Po = (io) / ( e o) — 0.6 ( Parge fc Pellat 1990| ). The mass density of each body is 1.5 g 
cm -3 . To set the initial mass of the annulus, M , we extend the Minimum Mass Solar 
Nebula to the Kuiper Belt and integrate the surface density distribution for solid particles, 
S = So(a/cto) across the 6 AU annulus. The dust mass is then M m , n « 0.25 So Me at 
32-38 AU for a® = 1 AU. Most Minimum Mass Solar Nebula models have So = 30-60 g 
cm" 2 QWeidenschilling 1977) ; |Hayashi 19811 ; |Bailey 19941 ), which sets M min w 7-15 M E . We 



consider models with M = 1-100 Me to allow for additional uncertainty in S . Table 1 
compares input parameters for Kuiper Belt models with initial conditions at 1 AU (see also 
W5gj ). 



Our success criteria are based on available observations of KBOs. The present day 
Kuiper Belt contains (a) at least one object (Pluto) with radius ^ 1000 km and (b) ~ 
70,000 objects with radii exceeding 50 km ( |Jewitt fc Luu 199~5| ; |Jewitt et al. 1996] , 1998). 



A successful KBO calculation must satisfy both observed properties in ^ 100 Myr. We 
quantify these criteria by defining r max as the radius of the largest object and r 5 as the 
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radius where the cumulative number of objects is No > 10 5 : a successful simulation has 
r max ^ 1000 km and r 5 J> 50 km at r ^ 100 Myr. To provide another characteristic of these 
models, we define r 95 % such that 95% of the mass is contained in objects with < r 95 %. In 
models with a long runaway growth phase, the largest objects contain most of the mass 
and Tq5% ~ r%. Models with a limited runaway growth phase leave most of the mass in 
small objects, so r 95 % <C r 5 . We end calculations with velocity evolution when r max exceeds 
~ 1000 km. To evaluate the dependence of runaway growth on fragmentation, we extend 
calculations without velocity evolution to 5000 Myr or to when r max exceeds ~ 2000-3000 
km. 



3.2. Models Without Velocity Evolution 

To isolate how fragmentation changes the growth of KBOs, we begin with constant 
velocity calculations of accretion. The initial size distribution has qo = 3, 5 = 1.4, and a 
maximum initial radius of r$ = 80 m. The initial velocities of hi = 4 m s _1 and Vi = 2.1 m 
s _1 correspond to an equilibrium model with eo = 10~ 3 ( [Hornung et al. 1985| ). We use the 



Davis et al. (1985) fragmentation algorithm and assume that the collision fragments receive 
a kinetic energy per unit mass equal to one-half of the square of the relative velocity of the 
colliding bodies, Vy. The bodies are strong, with a tensile strength S = 2 x 10 6 erg g _1 . 
Tables 1 and 4 summarize the initial conditions and results for models with Mq = 1-100 
M E , e = 10~ 3 , and the coefficients of restitution, (ci,c 2 ) = (10~ 5 ,10~ 5 ) and (10~ 2 ,10~ 3 ). 

Figure 3 shows how Nq evolves with time for M = 10 Me, eo = 10~ 3 , and {c\,C2) = 
(10 _5 ,10~ 5 ). The low coefficients of restitution eliminate rebound collisions. Roughly half 
of the most massive objects experiences at least one collision by r pa 16 Myr, when the 18 
largest bodies have m 1 km. These bodies reach r« ~ 10 km at r m 135 Myr and grow to 
100 km sizes at 255 Myr. The growth rate then increases considerably due to gravitational 
focusing. Runaway growth ensues. The largest planetesimals reach r max w 200 km at r pa 
265 Myr; r max exceeds 1000 km only 11 Myr later. At r pa 280 Myr, a single runaway body 
with r max pa 2000 km begins to sweep up lower mass planetesimals and contains nearly all 
of the mass in the annulus a few Myr later. 

The cumulative size distribution of small objects with <^ 100 m slowly approaches a 
power law with q = 2.25 throughout the evolution. The distribution is shallow, q < 2, for 
t <^ 20 Myr as growth produces many 50-500 m objects. The distribution steepens as the 
largest bodies grow past 1 km and settles at q — 2.25 for r ^ 150 Myr. This distribution is 
steeper than the theoretical limit of q = 2.5 ( Pohnanyi 1969 ; Williams fc Wcthcrill 1994]) , 



because growth and fragmentation never reach equilibrium. 
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Calculations with rebound collisions yield similar results. With larger coefficients of 
restitution, (01,02) = (10~ 2 ,10~ 3 ), rebound collisions occur for <^ 5 m. These collisions 
produce fragments but no mergers into larger bodies. As a result, the size distribution is 
initially very steep for r« <^ 5 m and very shallow for ^ 5 m (Figure 4). As the largest 
objects grow from 1 km to 10 km at r = 16-135 Myr, bodies with rj ^ 5 m are swept up 
by the large bodies and replaced by collision fragments. This process smooths out the size 
distribution at r\ ~ 1-20 m, although the slope still changes at r\ ~ 5 m. 

The evolution of the largest bodies is unaffected by rebounds. These objects reach sizes 
of 100 km only slightly later than models without rebounds, 258 Myr vs 255 Myr, and then 
begin to runaway from lower mass planetesimals. Several objects reach radii of 1000 km in 
another ~ 20 Myr and then begin to consume all of the material left in the annulus. 

Fragmentation is not important in any of these low eccentricity calculations. 
Catastrophic fragmentation produces no mass loss, because the velocities remain artificially 
low. These models lose only l%-3% of their initial mass due to cratering. The timescale to 
grow into a 1000 km object is nearly identical to models without fragmentation (see KL98): 
Pluto forms in r P ps 276-280 Myr (M /10 Me)" 1 for r = 80 m. Calculations with larger r 
should have smaller growth times (see KL98). 

3.3. Models With Limited Velocity Evolution 

To understand how fragmentation changes the velocity dispersions of colliding 
planetesimals, we consider calculations with 'limited' velocity evolution. As in Davis et 
al. (1985, 1994), we assume that collision fragments receive a fixed fraction f KE of the 
center-of-mass impact energy, but consider no other changes to the kinetic energy of the 
bodies. For most collisions, this assumption produces fragments with velocities larger than 
the relative collision velocity V^-. To conserve kinetic energy, the merged bodies have a lower 
velocity dispersion after the collision. This redistribution of velocity mimics dynamical 
friction and should lead to more rapid growth than models with no velocity evolution. Table 
5 summarizes results for eo = 10~ 3 , So = 2 x 10 6 erg g _1 , and coefficients of restitution, 
(ci, c 2 ) = (10~ 2 , 10~ 3 ), for our standard mass distribution with q = 3, 5 = 1.4, and r = 80 
m. 

Figure 5 shows how Nq evolves with time for Mo = 10 Me- During the first 10-20 
Myr, the size distribution has two power laws, q ~ 4 for <^ 5 m and q ~ 1.25 for ^ 5 m. 
This break in the slope is due to rebounds, which prevent growth of the smallest objects. 
As the largest objects grow to sizes of 1-10 km, collision debris adds to the population of 
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5 m objects. The slope of the size distribution is then q px 2.5 for <^ 100 m. For large 
objects, q approaches 2.5 once r max exceeds 100 km. Once r max ^ 1000 km, the complete 
number distribution follows a power law with q = 2.5, except for a small kink at ~ 1 km 
whose amplitude decreases with time. 

Limited velocity evolution eventually produces larger objects on shorter timescales 
than models without velocity evolution. In the first 20-40 Myr, collisions between the 
largest bodies, Tj pa 100-500 m, add material with high velocity dispersion to batches with 
Ti ~ 1-10 m. The horizontal velocity dispersions of these low mass objects increase from 4 
m s _1 to 10-30 m s _1 in 17 Myr and reach a roughly constant value of ~ 20 m s _1 after 
200 Myr (Figure 6). This large velocity dispersion initially slows down the growth of the 
largest objects relative to models with no velocity evolution. It takes this model ~ 1 Myr 
longer to produce objects with r\ ~ 1 km, but these objects then have velocities ~ 1 m 
s _1 smaller than the initial velocity of the planetesimal swarm (see Figure 6). This smaller 
velocity dispersion enhances the growth of the largest objects once gravitational focusing 
becomes important. The model produces 10 km objects in 120 Myr and 100 km objects in 
202 Myr. This rapid evolution leads to runaway growth and the production of a few 1000+ 
km objects after another 15 Myr. 

Figure 7 compares r max in models with no velocity evolution (thin solid lines) and 
models with limited velocity evolution (thick solid lines). The annuli have initial masses 
M = 1, 10, and 100 Me- These models are nearly indistinguishable for r max <^ 10 km. 
The limited velocity evolution models then produce larger objects on shorter timescales. 
These models do not have a mass-dependent velocity evolution such as viscous stirring or 
dynamical friction, so the timescale to produce 1000 km objects is still inversely proportional 
to the mass in the annulus, r P = 216 Myr (M /10 M^)" 1 for e = 10" 3 . 

Other than the accelerated growth rate at late times, there is little difference between 
models with limited velocity evolution and models with no velocity evolution (see Tables 
4-5). Both sets produce comparable numbers of large KBOs, r { ^ 100 km, that contain 
most of the mass in the annulus when the first Pluto forms. These models also lose 
~ l%-3% of their initial mass due to cratering. The mass loss depends solely on the 
total kinetic energy of all the bodies, which is a constant in these models. Catastrophic 
fragmentation still produces no mass loss, because we assume relatively strong objects with 
S = 2 x 10 6 erg g -1 . We describe models with weaker bodies in the next section. 
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3.4. Models with Velocity Evolution 



We now consider a complete coagulation calculation with fragmentation and velocity 
evolution. In addition to velocity redistribution from fragmentation, these models include 
(i) gas drag, (ii) dynamical friction and viscous stirring from long-range (elastic) collisions, 
and (iii) dynamical friction and viscous stirring from short-range (inelastic) collisions 
(see the Appendix of KL98). To approximate the disappearance of the solar nebula on 
reasonable timescales, we assume that the gas density decays exponentially with time, 
p g {t) = p gfi e-*/ T « where p gfi = 1.18 x 1(T 9 (a/lAU)~ 21 / 8 (M /4M B ) g cm" 3 (WS93 and 
references therein). The influence of gas drag on the planetesimals thus decreases with an 
e-folding time of r g , which we set at r g = 10 Myr ( [Hartmann et al. 1998| ). We adopt 5 = 
1.4, go — 3, and tq = 80 m. The initial velocities are hi = 4.0 (eo/10~ 3 ) m s _1 and Vi = 2.1 
(eo/10~ 3 ) m s -1 . Tables 6-7 summarize the initial conditions and results for models with 
M = 1-100 M E and e = 10~ 4 to 10~ 2 . 

3.4.I. A Standard Model with e = 10~ 3 and S = 2 x 10 6 erg g' 1 

Figure 8 shows how Nq and hi evolve for Mq = 10 Me, eo = 10~ 3 , and So = 2x 10 6 erg 
g _1 . Inelastic collisions act rapidly to circularize the orbits of objects with r« ^ 5 m, which 
decrease in velocity from hi = 4 m s _1 to hi = 1-2 m s _1 in 5 Myr. Larger bodies have 
less frequent collisions; their orbits circularize on longer timescales. Small bodies with <^ 
5 m collide frequently, but rebounds prevent these bodies from circularizing their orbits. 
Instead, collision fragments increase the velocity dispersion of the lowest mass bodies. 
These processes produce an inverted velocity distribution, where bodies with r; « 5 m have 
smaller velocities than objects with n <^ 5 m and ^ 10 m. 

The steady damping of particle velocities enhances growth of the larger bodies relative 
to calculations with limited velocity evolution. The largest objects grow slowly for 8.5 Myr, 
when there are ~ 500 objects with = 1 km. The size distribution then has three main 
features: (i) a pronounced fragmentation tail with a modest velocity dispersion, hi ~ 1-3 m 
s" 1 ; (ii) a transition region with pronounced kink in the size distribution; and (iii) a group 
of rapidly growing bodies with low velocities, hi = 0.02 m s _1 (Figure 8; left panel). These 
features remain prominent from 1 km to 100 km in only 3.5 Myr. 

During runaway growth, dynamical friction and viscous stirring begin to increase 
particle velocities (Figure 8; right panel). This phase ends when r max reaches ~ 300 km at 
15 Myr. The velocity dispersions of the small to intermediate mass bodies are then large 
enough, ~ 10 m s _1 , to reduce gravitational focusing factors considerably. The largest 
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objects enter a steady growth phase, where their radii grow slowly with time. There are 2 
Charon-sized objects with r, ^ 500 km at 19 Myr; 50 "Charons" at 25 Myr, and ~ 150 
Charons at 36.5 Myr when the first "Pluto" appears. The size distribution then closely 
follows a power law with q = 2.5 for <^ 30 m and a steeper power law with q = 3 for = 
1-1000 km. The modest velocity dispersion of the largest objects, hi ~ 0.5 m s -1 , maintains 
steady growth for another 60 Myr: 4 objects have r { = 1450 km at 50 Myr; 8 objects have 
Ti = 2000 km at 100 Myr. 

Figure 9 illustrates how the largest objects grow as a function of Mq. All models 
begin with slow growth, where inelastic collisions reduce the velocity dispersions of the 
intermediate mass objects. Slow growth lasts longer as M decreases, because the growth 
rate depends on the collision rate. Viscous stirring prolongs the slow growth phase of small 
Mq models by counteracting collisional damping earlier in the evolution. Slow growth thus 
produces more 1-2 km objects in models with smaller Mq. Runaway growth eventually 
turns 1 km bodies into 100 km objects. This phase ends when bodies with r\ <^ 1 km have 
velocity dispersions exceeding 10-20 m s -1 . The large gravitational focusing factors that 
began runaway growth are then smaller by several orders of magnitude. The growth then 
returns to a steady phase where the largest bodies gradually approach radii of 1000 km. 

In these models, cratering also acts to damp runaway growth. As noted in Sect. 
3.3, collision debris increases the velocities of the smallest objects. Viscous stirring and 
dynamical friction enhance this evolution as the largest objects grow beyond 1 km. The 
velocity dispersions of small objects increase from ~ 20 m s" 1 at the end of runaway growth 
to ~ 100 m s _1 at 100 Myr. The largest objects then have ps 2000 km. Collision debris 
produces ~ 30% of this increase; dynamical friction and viscous stirring are responsible for 
the rest. 

Although all models with cq = 10~ 3 make at least one Pluto, only annuli with Mq ~ 
10 Me meet both success criteria. In each calculation, r 5 steadily increases with time 
during the slow growth phase (Figure 10). The number of 50 km radius KBOs increases 
dramatically during runaway growth, when gravitational focusing factors are large. At the 
end of runaway growth, r 5 begins to decrease with time as the largest bodies try to separate 
themselves from the rest of the mass distribution. This evolution is modest and short-lived, 
because the runaway growth phase ends quickly. The number of 50 km radius KBOs then 
approaches a roughly constant value which grows with increasing Mq. Based on Figure 10, 
annuli with Mq <^ 3 Me produce too few 50 km radius KBOs compared with our success 
criterion; annuli with Mq ^ 30 Me probably produce too many. Annuli with 3 Me <^ Mq ^ 
30 Me can produce 1 or more Plutos and roughly 10 5 KBOs in £ 100 Myr. 

This 'standard model' has several important input parameters, including the initial 
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eccentricity, the initial size distribution, and the intrinsic strength of the bodies. To 
understand how the evolution depends on these parameters, we now consider variations 
on the standard model. We begin with a discussion of e and then describe models with 
different input strengths and size distributions. 

3.4.2. Models with e = 1(T 4 and 1(T 2 for S = 2 x 10 6 erg g' 1 

The velocity dispersion of the planetesimal swarm affects growth primarily through 
gravitational focusing. Large bodies grow faster when velocity dispersions are small and 
gravitational focusing factors are large. Despite small e , all of our calculations begin in 
the "high velocity regime," where gravitational focusing factors are near unity. Collisional 
damping and dynamical friction reduce velocity dispersions to the low velocity regime where 
growth is rapid. Planetesimals with small eo can reach this regime first, so we expect that 
smaller e leads to more rapid growth. 

Calculations with eo = 10" 4 and 10~ 2 confirm these general considerations. In models 
with eo = 10~ 2 , the larger initial velocities reduce collisional damping compared to models 
with e = 10~ 3 . Gravitational focusing factors thus increase slowly. Additional cratering 
debris also counteracts collisional damping and combines with dynamical friction and 
viscous stirring to keep velocity dispersions large. As a result, models with eo = 10~ 2 
experience a prolonged linear growth phase. This phase is ~ 10 times longer than for 
e = 1CT 3 models. Slow growth ends when the largest objects reach the sizes needed to 
produce modest gravitational focusing factors. A short runaway growth phase eventually 
produces several Plutos. The final size distribution has a larger fraction of mass in more 
massive objects, as indicated by a large value for rgs%, but has roughly comparable numbers 
of 50 km radius KBOs (Table 6). These models also lose more mass to dust, ~ 50%, 
compared to the l%-3% lost in models with e = 10~ 3 . 

Growth is rapid in models with eo = 10~ 4 , because collisional damping quickly reduces 
the velocity dispersion. This change does not have a dramatic effect on the growth of 
KBOs, because damping quickly reduces particle velocities to the low velocity limit. The 
damping time is then independent of velocity. Viscous stirring is also less effective. In our 
calculations, the growth times for e = 10~ 4 models are a factor of ~ 2-10 smaller than 
e = 10~ 3 models (Table 6). Models with e ^ a few x 10~ 5 probably have similar growth 
times, but we have not investigated this possibility in detail. 

Figure 11 compares the size distribution near the end of our calculations for three 
different values of e . All models produce two power law size distributions, with q = 2.5 for 
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small objects and q = 3 for large objects. These power laws are connected by a transition 
region which moves to larger radii with increasing e$. The fragmentation population thus 
extends to larger radii with larger e Q . In contrast, the merger population (the steep power 
law at larger radii) extends to smaller radii with smaller e$. This feature of the calculations 
can be tested directly with observations ( |Kenyon fc Luu 1999|) 

3.4.3. Models with Weaker Bodies 

The intrinsic strength So of a planetesimal changes the growth rate by setting the 

impact energy Q d needed to disrupt colliding bodies. In our disruption model, Qd depends 

on the sum of So and the gravitational binding energy of a colliding pair of planetesimals, 

4:irGpoRl/15 (eq. [All] and [A13]). The gravitational term is small compared to So when 

1/2 

the combined radius of the merged object is R c <^ 0.035*0 km. The maximum radius of our 
initial size distribution, r = 80 m, falls below this limit for So <^ 25 erg g _1 . For collisions 
between equal mass bodies, we can then derive a simple expression relating the minimum 
strength necessary for growth and the velocity dispersion, 

So,min « 6 x 10 4 erg g" 1 ( Vij V (1) 

V 100 m s 1 / 

At the start of our calculations, this result yields 

So^n-SOOergg- 1 ^^) 2 . (2) 

Bodies with eo ^ 0.1 initially grow for the standard case with 5*0 = 2 x 10 6 erg g -1 . As these 
objects grow to 1 km sizes, viscous stirring and dynamical friction increase the velocities 
of small objects. Once the velocities reach the threshold set by Equation (1), catastrophic 
disruption produces debris that is lost from the calculation. This process should limit 
the growth of the largest objects by reducing the reservoir of small bodies available for 
accretion. We thus expect the maximum size of KBOs to depend on So- 

To test these ideas in detail, we consider models with various S . Calculations with S 
= 10 to 2 x 10 6 erg g" 1 for e = 10" 4 and S = 10 4 to 2 x 10 6 erg g" 1 for e = 10" 3 allow 
velocity evolution to increase particle velocities up to the disruption threshold. Models with 
So ^ 2 x 10 6 erg g -1 and eo ^ 10~ 3 lose too much mass in the early stages to be of much 
practical value. Models with So ^ 2 x 10 6 erg g _1 are identical to models with So = 2 x 10 6 
erg g- 1 . 
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Table 7 summarizes our results at 50 Myr and 100 Myr for M = 10 Me and several 
different values of 5. For both values of eo, stronger objects can grow to larger sizes at 
100 Myr (see Figure 12). In each model, accretion and collisional damping lead to a short 
runaway growth phase that produces 100+ km objects with low velocities (hi <^ 0.1 m s _1 ), 
but leaves most of the initial mass in 0.1-1 km objects with much larger velocities (hi ~ 
3-10 m s _1 ). Dynamical friction and viscous stirring then increase the velocities of these 
small objects to the disruption threshold. The timescale to reach this threshold increases 
with So] r max also increases with So as indicated in Table 7 and Figure 12. When all small 
bodies have been disrupted, the maximum radius is nearly independent of 5: 



log r max « 2.45 + 0.22 log S (3) 

for 10 erg g _1 < So < 10 4 erg g _1 and eo = 10 -4 and 10 -3 . We did not run models with 
larger So to the disruption threshold. Future calculations will allow us to see whether r max 
reaches a threshold value at large So or continues to increase as indicated by Equation (3). 

Although r max depends on So, both r 5 and the slope of the final size distribution at 
large radius are independent of So- The radius limit for 10 5 KBOs has a small range, 
r§ fa 45-60 km, for any combination of 5 and So (Table 7). Figure 13 shows the evolution 
of the size distribution for So = 10 erg g _1 . As in Figure 8 and Figure 11, accretion 
produces a fragmentation tail with q = 2.5 at small radii and a steeper power law with 
q = 3 at large radii. The q = 3 power law persists throughout the catastrophic disruption 
phase; the fragmentation tail evolves into a very steep power law with q = 4. This 
behavior of the fragmentation tail occurs because larger objects initially experience more 
disruptive collisions than do smaller objects. Catastrophic disruption adds high velocity 
fragments from larger objects to the smaller mass bins; kinetic energy from this debris and 
viscous stirring gradually push smaller and smaller objects over the disruption threshold. 
Eventually, the smallest objects in our grid, n ~ 1 m, reach the disruption threshold and 
are slowly removed from the calculation. 



3.4-4- Models with Different Size Distributions 



The initial size distribution is one of the most uncertain input parameters of our 
coagulation models. The growth of 1-10 m or larger bodies from interstellar dust grains 
is poorly understood. Predicted sizes for conditions in the outer solar system range 
from ^ 1 m up to several hundred km depending on details of both microscopic and 
macroscopic physics (e.g., |Goldreich fc Ward 1973| ; Tremaine 1990| ; |Cuzzi et al. 1993a] , 
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1993b; [Bailey 1994| ; |Weidenschilling 19971 [Boss 19971 |Wurm fc Blum 1998|) . To test the 
sensitivity of our results to the initial conditions, we consider models with (a) q = 4.5, 
(b) q Q = 1.5, and (c) N c = Const 8(r — r ). These models produce final size distributions 
that are indistinguishable from our standard model with go — 3. Both r max and r§ are also 
independent of the initial size distribution. The time to runaway growth and the timescale 
to produce one Pluto increase with the amount of material initially in the smallest objects, 
because collisional damping is more effective when there is a large reservoir of small objects. 
We derive Tp = 30 Myr for q = 4.5, Tp = 37 Myr for q Q = 3, Tp = 42 Myr for q = 1.5, and 
t p = 49 Myr for N c = Const 5(r -80 m). 

To test the importance of the mass spacing factor, we recomputed these models for 5 
= 1.25. Our results for (a), (b), and (c) in the preceding paragraph are indistinguishable 
from results with 5 = 1.4, except that the timescale to produce Pluto decreases by ~ 5%. 
The model with 5 = 1.25 and our standard initial size distribution, q = 3, produced a 
single large body that ran away from the rest of the large bodies. The final radius of this 
single object is ~ 50% larger than r max for other calculations with M = lOM^. Several 
tests indicate that forming a single large body is a stochastic process sensitive to 5: small 5 
models produce such an object more often than large 5 models. We plan to investigate this 
possibility further in future studies. 



3.4-5. Other Input Parameters 

To conclude this section, we briefly comment on the sensitivity of our calculations 
to other input parameters listed in Table 1. The results described above are insensitive 
to factor of two variations in the particle mass density po, the relative gas velocity 77, 
the minimum velocity for cratering Vf, and the crushing energy Q c . We suspect that 
factor of ten variations in 77, Vf, and Q c will also have no impact on the results. Larger 
variations in p Q would probably change the variation of r max with Sq in Table 7, but we 
have not investigated this possibility. Increasing the fraction of kinetic energy imparted to 
fragmentation debris fxE decreases the time needed to produce 100+ km objects^. The 
growth time decreases by ~ 10% for fxE = 0.1 and ~ 20% for fxE = 0.2. Variations in the 
input j ke do not change the final size distribution or the number of KBOs as a function of 
Mq. This parameter thus has less impact on the results than either eo or So- 



2 In contrast, use of the WS93 fragmentation algorithm increases the time needed to produce 100+ km 
objects, although this increase is <; 10%-20% 
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3.5. Limitations of the Models 



We summarized the major uncertainties of our planetesimal calculations in KL98 
and will repeat important points for the current models here. Our choice of a single 
accumulation zone does not allow us to follow the evolution in semimajor axis of a 
planetesimal swarm flSpaute et al. 1991 ; |Weidenschilling et al. 1997|) . Although multi-zone 



calculations are important for understanding migration and other long-term aspects of 
planetesimal evolution, single-annulus models are a reasonable first approximation for 
the early evolution of KBOs. The growth of large nearby bodies, such as those that will 
merge to form Neptune, should modify the velocity evolution of KBOs once these bodies 
reach sizes much larger than 1000 km (sec [Morbidclli fc Valsccchi 1997| ). For most solar 
nebulae, these long-distance interactions should remain small until Pluto forms beyond 30 
AU (KL98; see also |Lazzaro et al. 1994| ; |Roques et al. 1994| ; |Morbidelli fc Valsecchi 1997] ; 



Ward fc Hahn 1998|). Calculations now underway will test this assertion in more detail. 



Our relatively coarse mass grid probably overestimates the timescale to produce KBOs 
and Pluto by ~ 5%-10% (KL98). The delay in runaway growth relative to a 5 = 1.25 
model is 3%-5% for 5 = 1.4 and 5%-10% for 5 = 2.0. Although these delays are small 
compared to the overall uncertainties in our algorithms, calculations with 5^2 rarely 
produce one or more isolated bodies that grow much faster than smaller bodies. The 
better mass resolution of S ^ 1.4 calculations allows more mass batches to satisfy the 
isolation condition, defined in the Appendix, which leads to more accurate calculations of 
the cumulative mass distribution and the velocity evolution of the lowest mass objects. 

The lack of a rigorous treatment of gas dynamics probably has little impact on our 
results. Gas drag removes ^ 1% of the initial mass from these models. Gas accretion by 
large bodies is also insignificant. The minimum radius needed to capture gas from the disk 
is n ~ 1000-2000 km for typical temperatures of 50-100 K at 30-50 AU ( |Beckwith eTaJ. 



1990|; psterloh fc Beckwith 199~4]). Our models reach this limit on timescales exceeding the 



lifetime of the gaseous disk, so we expect little gas accretion by Kuiper Belt bodies. 

Our calculations probably overestimate the amount of mass lost to dust. At both 1 
AU and 35 AU, losses from catastrophic fragmentation, cratering, and gas drag grow with 
increasing 6. In large 5 models, small delays in the formation of runaway bodies allow 
dynamical friction and viscous stirring extra time to increase the velocity dispersions, and 
hence mass loss, of the smallest objects. These effects are probably small, ~ 10%-20%, for 
most of our calculations. 

Our choice of the initial size distribution has little impact on our results. Models with 
go = 1-5 and go — 4.5 produce final size distributions very similar to those for calculations 



17 



with q = 3. The timescale to produce 1000 km objects lengthens as q decreases. The 
'equilibrium' size distribution with q = 3 is similar to the observed size distribution of 
interstellar dust grains, q ~ 2-3 (e.g, |Kim et al. 1994]; |Li fc Greenberg 1997| and references 



therein). Similar size distributions are derived from calculations of the growth of very small 
particles using measured sticking efficiencies (e.g., [Wurm fc Blum 1998|) . 

Aside from the timescale to produce 1000 km objects and the amount of mass lost to 
gas drag and fragmentation, the initial eccentricity distribution probably also has little 
impact on our conclusions. The number of 50 km radius KBOs is not sensitive to eo; ^95% 
increases with eo only for eo ^ 10~ 2 (Table 6). The slope of the final size distribution is also 
insensitive to eo- 

We suspect that initial eccentricities outside the range considered here are unrealistic. 
Viscous stirring and gas drag appear to set a lower limit on the eccentricity of small objects, 
e ~ 10 -5 . Models with eo <^ 10~ 4 should thus closely resemble eo = 10~ 4 models. For 
e ^ 10 -2 , collisions lead to substantial fragmentation and mass loss from the numerical 
grid. Circularization is then less effective and growth is very slow (see Table 6). These 
calculations probably poorly approximate reality: dynamical friction between 1-100 m 
objects and smaller dust particles not included in our grid should reduce the velocities of 
1-100 m objects on very short timescales^. We expect that these calculations would then 
more closely resemble models with smaller eo, if dust particles can grow back into 1 m 
bodies. 

The final limitation of our model is the fragmentation algorithm. We adopt an 
energy-scaling fragmentation law, because other types of models have not been developed 
and tested for the low velocity conditions appropriate in the Kuiper Belt. Other scaling 
models seem to yield better results for main belt asteroids than do energy-scaling models, 
but these models assume that the collision can be approximated as a point-like impact on a 
large body (Pavis et al. 1994 ; Durda et al. 1998 ; [Ryan fc Melosh 1998|) . This assumption 



is quite good for the high speed impacts, ^ 100-500 m s _1 , of strong, dense objects in the 
inner solar system. Point-like impacts are probably rare for the lower velocity collisions of 
weaker, low density objects like KBOs. Our consideration of KBOs with a large range of 
intrinsic strengths suggests that an improved fragmentation model would not change our 
results significantly. 



3 This problem does not occur in models with small eo, because the mass loss is <; 1% of the initial mass. 
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4. DISCUSSION AND SUMMARY 



In KL98 and this paper, we have developed a time-dependent planetesimal evolution 
code to calculate the formation of KBOs in a single annulus outside the orbit of Neptune. 
The computer program includes coagulation with realistic cross-sections, energy-scaling 
algorithms to treat cratering and disruption, and velocity evolution using the statistical 
expressions of Hornung et al. (1985). Our numerical solutions match standard analytic test 
cases and generally reproduce the results of other accretion and collision calculations (e.g., 
WS93; |Davis fc Farinella 1997Q . 



Our calculations demonstrate that plausible models can satisfy current observations of 
the Kuiper Belt. Several Plutos and ~ 10 5 50 km radius KBOs form in Minimum Mass 
Solar Nebulae with eo ~ 10~ 3 on timescales of 20-40 Myr. Growth is more rapid in more 
massive nebulae and in planetesimal swarms with smaller initial velocities. The formation 
time is less sensitive to the initial size distribution, the intrinsic strength of KBOs, and 
other input parameters listed in Table 1. 

Each Kuiper Belt model yields a cumulative size distribution with two main features. 
Objects with ^ 0.1 km follow Nq oc r~ q with q = 2.5, as expected for collision fragments 
( Dohnanyi 1969 ; Williams fc Wcthcrill 1994 ). Larger objects with ^ 1-10 km follow a 



q ~ 3 power law over several orders of magnitude in radius. These slopes do not depend on 
M , e , So, fxE-, and q , among other input parameters. Kenyon & Luu (1999) compare 
these results with observations. Here, we note that the q m 3 power law for large bodies is 
identical to the observed slope, 3 ± 0.5 QJewitt et al. 1998] ; [Luu fc Jewitt 19981 ; but see also 
Gladman et al. 1998] ). 



Fragmentation and velocity evolution are important components in the formation of 
present day KBOs. Fragmentation produces a large reservoir of small bodies that damp the 
velocity dispersions of the large objects through dynamical friction. These processes allow 
a short runaway growth phase where 1 km objects grow into 100 km objects. Continued 
fragmentation and velocity evolution damp runaway growth by increasing the velocity 
dispersions of small objects. This evolution leaves ~ l%-2% of the initial mass in 50 km 
radius KBOs. The remaining mass is in 0.1-10 km radius objects. Fragmentation will 
gradually erode these smaller objects into dust grains that are removed from the Kuiper 
Belt on short timescales, ~ 10 7 yr (see [Backman fc Paresce 1993] ; |Backman et al. 1995| ). 



Thus, 50 km radius KBOs comprise a small fraction of the original Kuiper Belt. 

Fragmentation also limits the size of the largest object in the Kuiper Belt. The 
maximum radius ranges from r max « 450 km for Sq = 10 erg g" 1 to r max ^ 3000 km for 
So = 2 x 10 6 erg g . Pluto formation sets a lower limit on the tensile strength of KBOs, 
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S > 300 erg g" 1 . 



These results suggest a refinement of our picture for KBO formation in the outer solar 
system. In KL98, we speculated that velocity perturbations due to the growth and outward 
migration of Neptune would limit the growth of KBOs at radii <^ 1000 km. Although this 
hypothesis is plausible (see, for example, [Vlalhotra 1993| ; |Morbidelli fc Valsecchi 19971) , our 
current models demonstrate that 50-1000 km radius KBOs form naturally at ~ 35 AU on 
timescales, Tp ~ 10-100 Myr, comparable to the Neptune formation time. Although a few 
objects can reach 2000+ km radii on timescales of 2-3 Tp, nearly all of the initial mass 
beyond 30 AU remains in small, 1 km radius objects that can be depleted by collisional 
disruption ( |Davis & Farinella 1997| ) or gravitational sculpting ( jHolman fc Wisdom 1993] ) or 
both on timescales exceeding 100 Myr. This evolution can account for the observation of 
50+ km KBOs in a currently small mass Kuiper Belt without intervention by Neptune. 

Finally, our new results further support the notion that KBOs will form in other 
solar systems. The dusty circumstellar disks detected in many pre-main-sequence stars 
suggest masses of 1-100 M E at distances of 30-100 AU (e.g., |Sargent fc Beckwith 1993| ; 
Beckwith fc Sargent 1996] ; see also |Close et al. 1997| ; [Hogerheijde et at. 1997| |hay et al. 



TW?] |Akeson et al. 1998|; [Stapelieldt et al. 1998]). KBOs with 50+ km radii can grow 



in this material as the central stars contract to the main-sequence if the disks are not too 
turbulent (see |Cuzzi et al. 1993a] , 1993b). Smaller, 1-10 km radius KBOs probably form in 
less massive pre-main-sequence disks. Our results also indicate that the growth of 100-1000 
km radius KBOs is accompanied by substantial dust production, ~ 0.1-1 Me, in models 
with M ~ 10-100 Me- This dust could be responsible for the ringlike structures observed 
in pre-main-sequence stars such as GG Tau ( |Roddier et al. 1996| ) and older main sequence 
stars such as e Eri (|Greaves et al. 1998| ) and HR 4796 ( [Jayawardhana et al. 1998] ; |Koerner 



et al. 19§8]). The less massive disks in a Lyr, a PsA, and (3 Pic may also harbor KBOs if 



the dust masses are reasonably close to the 'maximum' masses inferred for these systems 
( [Backman fc Paresce 1993] ). In all of these stars, the dynamics and mass distribution of 
dust may well provide useful constraints on the properties of presumed KBOs. We hope to 
explore this possibility in future studies. 
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Science. Comments from A. Cameron, F. Franklin, M. Geller, M. Holman, S. Starrfield, 
and J. Wood greatly improved our presentation. We acknowledge G. Stewart for clarifying 
details of the WS93 calculations. 
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A. APPENDIX 



A.l. Overview 

As described in KL98, our accretion model assumes that planetesimals are a statistical 
ensemble of masses in a cylindrical annulus of width Aa and height H centered at a radius 
a from the Sun. The particles have horizontal hi(t) and vertical Vi(t) velocity dispersions 
relative to an orbit with mean Keplerian velocity Vr (see [Lissauer fc Stewart 1993| ). We 



approximate the continuous distribution of particle masses with discrete batches having an 
integral number of particles n.j(i) and total masses Mj(t) ( |WS93| ). The average mass of a 
batch, TOj(t) = Mi(t)/rii(t), evolves with time as collisions add and remove bodies from the 
batch. 

To compute the evolution of particle numbers and velocities in KL98, we solved the 
coagulation equation and a set of velocity evolution equations for all mass bins k during 
a time step St. We assumed that bodies merge but do not fragment during collisions. 
We conserved kinetic energy in each collision and adopted a kinetic approximation to 
calculate velocity changes due to gas drag, dynamical friction, and viscous stirring. Our 
explicit algorithm for solving these equations reproduced standard tests and other published 
calculations. 

In this paper, we consider collisions that produce mergers and debris. The coagulation 
equation is then: 

5n k = 5t [eijAijUiUj - n k A ik n^ + 8n k j - Sn k , g d (Al) 



5M k = 8t < ijAjjiiiii ,in k - n k A ik nim k } + m k 5n k j - m k 5n Kgd (A2) 

where is the cross-section, = 1/2 for i — j and = 1 for i ^ j. The four terms in 
A1-A2 represent (i) mergers of rrii and rrij into a body of mass m k = rrii + rrij — vn e ^^ (ii) 
loss of m k through mergers with other bodies, (iii) addition of m k from debris produced 
by the collisions of other bodies, and (iv) loss of m k by gas drag. We consider below the 
mass lost to small fragments vn e ^j. The second term in A1-A2 includes the possibility 
that a collision can produce debris but no merger (rebounds). As in KL98, we calculate 
the "gravitational range" of the largest bodies - R g ^ = KiaRH,u mid + 2aej (|WS93|) - where 



Ki = 2 v^3 and Rn,ij = [(mi + rrij) /?>~M. & ] 1 ^ is the mutual Hill radius. As in WS93, the 
isolated bodies are the iV largest bodies that satisfy the summation, Y^r ax n iRg,i > Aa. 
These isolated bodies cannot collide with one another but can collide with other lower mass 
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bodies. 

As in KL98, we solve the complete set of evolution equations, A1-A2 above and A7-A8 
from KL98, using an explicit method that limits the time step automatically to prevent 
large changes in the dynamical variables. Section 2 of the main text compares calculations 
at 1 AU with results from WS93. In the rest of this Appendix, we describe fragmentation 
algorithms (A. 2) and updates to our velocity evolution algorithm (A. 3). 



A. 2. Fragmentation 

Algorithms for collision outcomes rely on comparisons between the kinetic energy of the 
two colliding planetesimals Q c ^ and the binding energy of the merged planetesimal Qb,ij- 
The binding energy usually consists of an intrinsic tensile strength and the gravitational 
binding energy, Q^j — Sq + Q g ,ij- This 'energy scaling' approach is sometimes replaced 
by other scaling laws to model the structure of the colliding bodies more accurately. For 
example, Housen and collaborators (see [Housen fc Holsapple 1990| ; |Housen et al. 1991| ; 



|Holsapple 1993| , 1994) describe a strain-rate scaling model to express the collision energy 
needed to disrupt a body in terms of its size and impact velocity. Davis et al. (1994) showed 
that both types of model can match observations of the asteroid belt in our solar system 
(see also Farinella et al. 1982j ; Pavis et al. 1989| ; |Williams fc Wetherill 1994] |Marzari et al 



[19971 ). 



In this paper, we adopt an energy scaling model. Energy scaling has two main 
advantages for calculating the fragmentation of Kuiper Belt objects. The wide use of this 
approach allows us to compare results with many previous calculations. Other models also 
seem inappropriate for the low velocity collisions anticipated in the Kuiper Belt. These 
models have been scaled for conditions in the inner Solar System, where collision velocities 
exceed 100 m s" 1 (see Pavis et al. 1994]) . The approximations made for these collision 
velocities fail when the impact velocity is smaller than 100-1000 m s -1 (e.g., [Housen 15 



[Holsapple 1990| , p. 239). We expect velocities of ^ 10 m s -1 for most collisions in the 
Kuiper Belt (KL98), and thus cannot apply the Housen et al. (1990) and other sophisticated 
models in our calculations. 

The output of any fragmentation algorithm is the mass of the merger remnant m& and 
the mass distribution of the fragments. We assume rrik = rrii + rrij — m ey ij, where m e ^j is 
the mass of fragments ejected from the merged planetesimal. We consider two approaches 
to compute m e ^. In the first case, we follow WS93 and set the impact velocity as 
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Vj, V = K + K% (A3) 

where is the relative velocity of the two colliding bodies (equation A12 of KL98) and 
= IGivrii + rrij)/(ri + Tj) is the mutual escape velocity. The center of mass fragmentation 
energy Q f ti j and the gravitational binding energy Q g ^j per unit mass are 



Qui = (^/2) rrii rnj V/^/K + mj) 2 (A4) 

and 

Q gtij = 0.6 G (mi + rrij)/R c (A5) 

where R c is the spherical radius of the combined body with mass mj + rrij and K2 = 0.5 is 
a constant. We further define the fragmentation energy Ef^j = (rrii + m j)Qf,ij- 

WS93 and most other studies assume that collisions produce debris when (i) the 
impact velocity exceeds a threshold velocity, Vj^j > Vf, and (ii) the amount of ejected 
mass exceeds a threshold value, m ei ij > 10 _8 (mj + rrij). The mass ejected in the collision 
is derived from simple energy considerations. The collision disrupts the merged object if 
the collision energy exceeds the binding energy, Qf t ij > Q g ,ij + So, where So is the impact 
strength. In most of their calculations, WS93 further require that the mass fragmented by 
the impact exceed half of the mass of the merged object, rrif^j > 0.5 (rrii + rrij), where the 
fragmented mass is ( |Greenberg et at. 1978| , 84): 



m f,ij = Qf,ij/Qc , (A6) 

and Q c is the crushing energy. If both conditions for "catastrophic disruption" are met, the 



collision produces a large fragment with mass (Fujiwara et al. 1977; see also [Fujiwara 1980 



Housen fc Holsapplc 1990|) : 



m L ,ij = 0.5 (mi + mj) (Q m ,ij/(Q g ,ij + Qb)) ' , (A7) 

and numerous smaller fragments. In all cases, m^^j < 0.5 (m 8 + mj). 

If the conditions for catastrophic disruption are not met, the collision produces a large 
body with m^ pa mi + rrij and numerous small fragments with total mass rrif^j (A6). WS93 
assume that the amount of the fragmented mass that escapes the merged body is 
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m 



where K 3 = 3 x 10 (cm s 
1978| ). The mass of the largest fragment in this case is 



K*m M V-£* (A8 
x ) 2 - 25 is a constant (see also |Gault et al. 1963] ; preenberg et a l. 



m L)ij = 0.2 m e ,ij 



(A9) 



In WS93, the mass distribution of the debris depends on the method of fragmentation. We 
simplify this procedure and adopt a single power law distribution for both methods: 



5n k = C k {m x h - m 2 b ) 



where b = (1 + '/'/.,,•'//', ,;) 1 , C k 



m 



mi 



{m k m k _i 



) 1/2 , and 



m 2 



{m k m k+1 



(A10) 



,1/2 



( preenberg et al. 1978|) . Tests show that this distribution produces results similar to those 
of WS93 when we assume that the velocity of the escaping fragments equals the relative 
velocity of the colliding planetesimals, (see section 2). 

As a comparison to the WS93 fragmentation algorithm, we consider the Davis et al. 
(1985; |Greenberg et al. 1978| , 1984; Pavis et al. 1994]) approach. We adopt their energy 
scaling formula and write the strength (in erg g _1 ) of a planetesimal as: 



S = Sn 



A-kK a G Po R 
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(All) 



where K4 is a constant. We adopt K4 = 1 in the absence of any useful estimates for Kuiper 
Belt objects; Davis et al. (1985, 1994) consider K A = 1-100 for collisional evolution in 
the asteroid belt (see also [Housen &; Holsapple 199U ). Davis et al. assume that the ejecta 
receive a fixed fraction fKEQf,ij of the impact energy and have a power law mass- velocity 
distribution, 



/(> v) oc (v/v c 



-ay 



(A12) 



where v c is a reference velocity. The impact energy needed to disrupt the planetesimal and 
accelerate 50% of the fragments to escape velocity follows from integration of (A12) and 
energy conservation ( |Davis et al. 1985] ): 



24 



V jKE{cev - 2) J 

Davis et al. assume that collisions produce debris when the impact velocity exceeds 
a threshold, Vj^j > Vf. We also adopt the WS93 threshold for the ejected mass, 
m e> ij > 10~ 8 (mj + rrij). The collision disrupts the colliding bodies if Q m ,ij > Qd- If this 
condition is met, the total mass lost is: 

m eM = 0.5 (mi + m,) (Q^ij/Qdf 5 ^ , (A14) 
while the mass of the largest fragment is the smaller of 

_ f 0.5 K + mj ) (Qm.v/Qd) 1 -** , A1 ^ 

mL ^ - \ 0.2 m^j (A15) 

If the collision does not disrupt the merged planetesimal, the mass in fragments is simply 
rrifjj from (A6). To estimate the fraction of this mass that escapes, we follow the Davis et 
al. (1985) derivation of Qd and integrate the mass distribution over velocity: 

m e,ij = 772 ( A16 ) 

The mass of the largest fragment is then m^^ = 0.2 as before. We calculate the mass 
distribution of the fragments as in (A10) and assume that all fragments have the same 
kinetic energy per unit mass. 

Unlike WS93, Davis et al. allow for 'rebounds,' where the colliding bodies do not merge 
into a single body. We follow Barge & Pellat (1992) and define the rebound velocity as 

V 2 2 (! ~ Cr) V 2 ( A)7) 

v reb,ij o y e,ij K^^') 

C R 

where cr is the coefficient of restitution. In most applications, cr takes on separate values 
for collisions with impact velocities above and below the fragmentation threshold: 

[ c 2 V Itij > V f 
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With these definitions, collisions with Vj^j < V re b ti j produce mergers and debris; collisions 
with Vi ti j > V re b t ij produce debris but no merger. We assume that rebounds conserve energy 
and that the velocity difference between bodies after the collision is 



{Vi - V^after = C R (Vi - Vj) before . (A19) 

Tests of the coagulation code indicate that this prescription yields results almost identical to 
a prescription where we assume that the particle velocities are not changed by the rebound. 

In addition to the tests described in section 2, we repeat the collisional evolution of 
KBOs described in Davis & Farinella (1997). With suitable modifications to our algorithms 
for the collision cross-section and velocity evolution, we reproduce their results to ~ 
10%-20%. 



A. 3. Velocity evolution 

Our procedure for treating the evolution of the horizontal hi and vertical v $ components 
of the velocity dispersion follows the kinetic formalism developed by Hornung et al. (1985) 
and Barge & Pellat (1990, 1991, 1992). In KL98, we adopted the WS93 expressions for 
dynamical friction and viscous stirring due to long range encounters and reformulated the 
Barge & Pellat (1990, 1991) expressions for inelastic collisions in terms of e and i (see 
equations A19-A22 of KL98). We verified that our algorithm for computing the velocity 
evolution achieved the equilibrium value for the ratio of the vertical to the horizontal 
velocity j3 = i/e = 0.6 ( Hornung et al. 1985|) and generally reproduced the velocity evolution 



of the WS93 Earth calculation. 

In this paper, we make two modifications to our treatment of the velocity evolution 
and clarify our treatment of velocity evolution at low velocities. We assume that rebound 
collisions do not contribute to dynamical friction and viscous stirring for inelastic collisions 
(equations A21 and A22 of KL98). With this assumption, the orbits of low mass objects do 
not circularize when their velocity dispersion exceeds the rebound velocity (equation A17 
above). This behavior slows down the growth time of massive objects, because the collision 
cross-section remains in the high velocity limit. 

Although the kinetic prescription of Hornung et al. (1985) accurately represents 
dynamical evolution at high particle velocities, it fails as the particle velocity approaches 
the Hill velocity v H>ij = VtR H ,ij Qlda 1990t |Barge fc Pellat 1991| ; WS93). Ida (1990) used 



A^-body calculations to derive the appropriate behavior at low particle velocities, V« < 2-5 
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v H,ij- We generally reproduce Ida's (1990) results if we adopt 




'2 2 \ ^ / J °2 1 

e i + e j \ ( di> VSi i 



e ij,lv J V dt j } r 



(A20) 



for viscous stirring and 



'Hp 2 \ f p 2 -4- p 2 \ /Hp 2 \ /Hi 2 \ / p 2 4- p 2 \ /Hi 2 ' 



dt J lv \ e%j,iv J \ dt J hv \ dt J [v \ eijjv J \ dt J hv 



(A21) 



for dynamical friction (see also Parge fc Pellat 1991| , |WS93| ). In both equations, the 'hv 



subscript indicates the appropriate high velocity expression from KL98 (equations A19 and 
A20) and is the value of e 2 + e 2 at Vij = Vi v vu,ij- These equations yield constant 
timescales for dynamical friction and the inclination component of viscous stirring ( |lda| 
1990|; Figure 10). The timescale for the eccentricity component of viscous stirring varies as 



e 2 at low velocity (|Ida 199CQ . 

We describe the sensitivity of the velocity evolution to the choice of V\ v in Sec. 2. Any 
Viv ^ 2 slows down velocity evolution significantly, because the most massive particles in 
our test calculations have ^ 0.3-0.5 vn,ij- The velocity changes are then ~ 25 times 
smaller for viscous stirring and ~ 500 times smaller for dynamical friction (see also [Ida 



T9"9"0| |Barge & Fellat 1991| , |WS93| ). Our results are not very sensitive to V iv for 2 £ V iv ^ 5: 



the velocity evolution slows down gradually as Vi v increases from 2 to 5 (see Sec. 2). We 
adopt the middle of Ida's (1990) range, Vi v = 3.5, in our Kuiper Belt calculations. For 
comparison, WS93 adopt Vi v = 2; Parge fc Pellat 1991| adopt V\ v = 0.8 for viscous stirring 
and V\ v = 2.6 for dynamical friction. 

In many of our Kuiper Belt calculations, the velocity of the largest bodies approaches 
zero due to collisional damping, dynamical friction, and fragmentation. To avoid divergences 
in the stirring rates from long-range encounters, we adopt a minimum horizontal velocity 
of 10~ 3 m s -1 and a minimum vertical velocity of 5.3 x 10 -4 m s -1 . This 'floor' to the 
velocity evolution maintains the appropriate for collisions with > V[ v vn,ij and allows 
larger time steps for the brief interval when growth is rapid and the viscous stirring rate is 
small. Tests show that the growth rate of the largest bodies does not depend on the values 
for these limits. Lower values for the 'floor' result in somewhat higher velocities for the 
small mass objects, but these changes are small. Once viscous stirring begins to dominate 
collisional cooling, the velocities of all bodies increase above our 'floor' values. 
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Finally, to reduce the computation time of our velocity evolution algorithm, we 
approximate the J, J, and K integrals in equations A19-A22 of KL98 with polynomial 
expressions derived using Mathematicc^\ Our results match the integral expressions to better 
than 1 part in 10 3 for /%■ = 0-1. Expressions ready for use in a computer calculation appear 
below. 

I r = 14.1439 + 0.0675783//% + /% * (3.86566 + /% * (-6.5623 

+ p, tj * (12.8084 + pij * (-9.93223 + 2.94993 * /%)))) (A21) 

I e = 2.32693 + 0.0337724//% + /% * (1.80115 + /% * (-3.3792 

+ p, tj * (5.78079 + pij * (-4.43822 + 1.32431 * /%)))) (A22) 

I z = -0.00137043 + /%* (0.121347 + /% * (7.27632 + (3 i3 * (6.32281 

+ /% * (4.98011 - 1.36123 * /%)))) (A23) 

J r . = -4.21891 + p, tj * (-0.762162 + /%• * (37.2078 + /%• * (-97.5891 
+ fa * (120.258 + pij * (-72.0619 + /% * (14.8515 
+ 1.57744* /%)))))) (A24) 

J e = -1.90771 + fa * (13.8618 + /%• * (-30.8032 + * (47.3479 
+ /% * (-49.4466 + /% * (32.4446 + /%• * (-11.7370 
+ 1.72263 * /%)))))) (A25) 

J z = 6.12664 + Pij * (-13.1003 + /% * (-6.3946 + /% * (50.1757 
+ fa * (-70.6068 + /%■ * (39.2953 + /% * (-2.86619 
- 3.37452 * /%)))))) (A26) 



4 Mathematica «3.0.2 ©1988-1997 Wolfram Research, Inc. 
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K r = 3.16452 + 0^ * (0.0149191 + Aj * (-1.73044 + Aj * (1.55516 
+ fa * (1.08438 + fa * (-3.55476 + Aj * (2.91763 
- 0.841079 * /%)))))) (A27) 



K Q = 2.96882 + 0.135176//% + A; * (7.13193 + Aj * (-24.1086 

+ ptj * (39.3287 + Pij * (-35.5738 + Aj * (15.0527 
+ p^ * (-0.327693 - 1.25529* Aj)))))) (A28) 



K z = 6.13154 + 0.13519/ P^ + pij * (-6.70136 + Aj * (4.49227 

+ pij * (-6.53516 + Pij * (19.317 + Aj * (-29.5368 
+ p^ * (21.0967 - 5.78973* A;)))))) (A29) 



29 



REFERENCES 

Akeson, R. L., Koerner, D. W., & Jensen, E. L. N. 1998, ApJ, 505, 358 

Bailey, M. 1994, In Asteroids, Comets, Meteors 1993, eds. A. Milani, M. DiMartino, and A. 
Cellino (Kluwer Academic Publishers, Dordrecht), pp. 443 - 459. 

Backman, D. E., & Paresce, F. 1993, in Protostars and Planets III, eds. E. H. Levy & J. I. 
Lunine, Tucson, Univ of Arizona, p. 1253 

Backman, D. E., Dasgupta, A., & Stencel, R. E. 1995, ApJ, 450, L35 

Barge, P., & Pellat, R. 1990, Icarus, 85, 481 

Barge, P., & Pellat, R. 1991, Icarus, 93, 270 

Barge, P., & Pellat, R. 1993, Icarus, 104, 79 

Beckwith, S. V. W., & Sargent, A. I. 1996, Nature, 383, 139 

Beckwith, S. V. W., Sargent, A. L, Chini, R., & Giisten, R. 1990, AJ, 99, 924 

Boss, A. P. 1997, Science, 276, 1836 

Close, L. M., Roddier, F., Hora, J. L., Graves, J. E., Northcott, M., Roddier, C, Hoffman, 
W. F., Dayal, A., Fazio, G. G., & Deutsch, L. K. 1997, ApJ, 489, 210 

Cuzzi, J. N., Dobrovolskis, A. R., & Champney, J. M. 1993, Icarus, 100, 102 

Cuzzi, J. N., Dobrovolskis, A. R., & Champney, J. M. 1993, Icarus, 106, 102 

Davis, D. R., Chapman, C. R., Weidenschilling, S. J., & Greenberg, R. 1985, Icarus, 62, 30 

Davis, D. R., & Farinella, P. 1997. Icarus, 125, 50 

Davis, D. R., Ryan, E. V., & Farinella, P. 1994, Planet. Space Sci., 42, 599 

Davis, D. R., Weidenschilling, S. J., Farinella, P., Paolicchi, P., & Binzel, R. P. 1989, in 
Asteroids II, edited by R. P. Binzel, T. Gehrels, & M. S. Matthews, Tucson, Univ. 
of Arizona Press, p. 805 

Dohnanyi, J. W. 1969, J. Geophys. Res., 74, 2531 

Duncan, M. J., Levison, H. F., & Budd, S. M. 1995, AJ, 110, 3073 

Durda, D. D., Greenberg, R., & Jedicke, P. 1998, Icarus, 135, 431 

30 



Farinella, P., Paolicchi, P., k Zappala, V. 1992, A&A, 253, 604 
Fernandez, J. A., & Ip, W.-H. 1981, Icarus, 47, 470 
Fernandez, J. A., & Ip, W.-H. 1984, Icarus, 58, 109 
Fujiwara, A. 1980, Icarus, 41, 356 

Gault, D. E., Shoemaker, E. M., & Moore, H. J. 1963, NASA TN, 1767 
Gladman, B., & Kavelaars, J. J. 1997, A&A, 317, L35 

Gladman, B., Kavelaars, J. J., Nicholson, P. D., Loredo, T. J., & Burns, J. A. 1998, AJ, 
116, 2042 

Goldreich, P., & Ward, W. R. 1973, ApJ, 183, 1051 
Greaves, J. S. et al. 1998, ApJ, 506, L133 

Greenberg, R., Wacker, J. F., Hartmann, W. K., & Chapman, C. R. 1978, Icarus, 35, 1 
Greenberg, R., Weidenschilling, S. J., Chapman, C. R., & Davis, D. R. 1984, Icarus, 59, 87 
Hartmann, L., Calvet, N., Gullbring, E., & D'Alessio, P. 1998, ApJ, 495, 385 
Hayashi, C. 1981, Prog Theor Phys Suppl, 70, 35 

Hogerheijde, M. R., van Langevelde, H. J., Mundy, L. G., Blake, G. A., & van Dishoeck, E. 
F. 1997, ApJ, 490, L99 

Holman, M. J., & Wisdom, J. 1993, AJ, 105, 1987 

Holsapple, K. A. 1993, Ann. Rev. Earth Planet. Sci., 21, 333 

Holsapple, K. A. 1994, Planet. Space Sci., 42, 1067 

Housen, K., & Holsapple, K. 1990, Icarus, 84, 226 

Housen, K., Schmidt, R. M., & Holsapple, K. 1991, Icarus, 94, 180 

Hornung, P., Pellat, R., & Barge, P. 1985, Icarus, 64, 295 

Ida, S. 1990, Icarus, 88, 129 

Ip, W.-H. 1989, Icarus, 80, 167 

Irwin, M., Tremaine, S., & Zytkow, A. N. 1995, AJ., 110, 3082 

31 



Jayawardhana, R. et al. 1998, ApJ, 503, L79 

Jedicke, P., k Metcalfe, T. S. 1998, Icarus, 131, 245 

Jewitt, D., k Luu, J. 1993, Nature, 362, 730 

Jewitt, D., k Luu, J. 1995, AJ, 109, 1867 

Jewitt, D., Luu, J., k Chen, J. 1996, AJ, 112, 1225 

Jewitt, D., Luu, J. X., k Trujillo, C. 1998, AJ, 115, 2125 

Kenyon, S. J., k Luu, J. X. 1998, AJ, 115, 2136 (KL98) 

Kenyon, S. J., k Luu, J. X. 1999, ApJL, submitted 

Kim, S.-H., Martin, P. G., k Hendry, P. D. 1994, ApJ, 422, 164 

Koerner, D. W., Ressler, M. E., Werner, M. W., k Backman, D. E. 1998, ApJ, 503, L83 

Kolvoord, R. A. k Greenberg, R. 1992, Icarus, 98, 2 

Lay, O. P., Carlstrom, J. E., k Hills, R. E. 1997, ApJ, 489, L917 

Lazzaro, D., Sicardy, B., Roques, F., k Greenberg, R. 1994, Icarus, 108, 59 

Levison, H. F., k Duncan, M. J. 1993, ApJ, 406, L35 

Li, A., k Greenberg, J. M. 1997, A&A, 323, 566 

Lissauer, J. J., k Stewart, G. R. 1993, In Protostars and Planets III, edited by E. H. Levy 
and J. I. Lunine, U. of Arizona Press, Tucson, 1061 

Lissauer, J. J., Pollack, J. B., Wetherill, G. W., k Stevenson, D. J. 1996. "Formation of the 
Neptune System." In Neptune and Triton, eds. D. P. Cruikshank, M. S. Matthews, 
and A. M. Schumann (U. of Arizona Press, Tucson, pp. 37 - 108. 

Luu, J. X., k Jewitt, D. 1998, ApJ, 502, L91 

Luu, J. X., Marsden, B., Jewitt, D., Trujillo, C. A., Hergenother, C. W., Chen, J., k Offutt, 
W. B. 1997, Nature, 387, 573 

Malhotra, R. 1993, Nature, 365, 819 

Malhotra, R. 1995, A J, 110, 420 



32 



Marzari, F., Farinella, P., Davis, D. R., Scholl, H., Bagatin, A. C. 1997, Icarus, 125, 39 

Morbidelli, A., & Valsecchi, G. B. 1997, Icarus, 128, 464 

Ohtsuki, K., & Nakagawa, Y. 1988, Prog Theor Phys (Suppl), 96, 239 

Ohtsuki, K., Nakagawa, Y., & Nakazawa, K. 1990, Icarus, 83, 205 

Osterloh, M. & Beckwith, S. V. W. 1994, ApJ, 439, 288 

Pollack, J. B. 1984, ARA&A, 22, 389 

Pollack, J. B., Hubickyj, O., Bodenheimer, P., Lissauer, J. J., Podolak, M., & Greenzweig, 
Y. 1996, Icarus, 124, 62 

Roddier, C, Roddier, F., Northcott, M. J., Graves, J. E. & Jim, K. 1996, ApJ, 463, 326 

Roques, F., Scholl, H., Sicardy, B., & Smith, B. A. 1994, Icarus, 108, 37 

Russell, S. S., Srinivasan, G., Huss, G. R., Wasserburg, G. J., & Macpherson, G. J. 1996, 
Science, 273, 757 

Ryan, E. V., & Melosh, H. J. 1998, Icarus, 133, 1 

Safronov, V. S. 1969, Evolution of the Protoplanetary Cloud and Formation of the Earth 
and Planets, Nauka, Moscow [Translation 1972, NASA TT F-677] 

Sargent, A. L, & Beckwith, S. V. W. 1993, Phy Tod, 46, 22 

Spaute, D., Weidenschilling, S. J., Davis, D. R., & Marzari, F. 1991, Icarus, 92, 147 

Stapelfeldt, K. R., Krist, J. E., Menard, F., Bouvier, J., Padgett, D. K., & Burrows, C. J. 
1998, ApJ, 502, L65 

Stern, S. A. 1995, AJ, 110, 856 

Stern, S. A. 1996, AJ, 112, 1203 

Stern, S. A., & Colwell, J. E. 1997a, AJ, 114, 841 

Stern, S. A., & Colwell, J. E. 1997b, ApJ, 490, 879 

Tremaine, S. 1990, In Baryonic Dark Matter, edited by D. Lynden-Bell & G. Gilmore, 
Kluwer, Dordrecht, p. 37 

Ward, W. R., & Hahn, J. M. 1998, AJ, 116, 489 

33 



Weidenschilling, S. J. 1977, Astrophys Sp Sci, 51, 153 
Weidenschilling, S. J. 1997, Icarus, 127, 290 

Weidenschilling, S. J., Spaute, D., Davis, D. R., Marzari, F., & Ohtsuki, K. 1997, Icarus, 
128, 429 

Wetherill, G. W. 1990, Icarus, 88, 336 

Wetherill, G. W., & Stewart, G. R. 1989. Icarus 77, 300 - 357. 
Wetherill, G. W., & Stewart, G. R. 1993, Icarus, 106, 190 (WS93) 
Williams, D. R., & Wetherill, G. W. 1994, Icarus, 107, 117 

Williams, I. R, O'Ceallaigh, D. P., Fitzsimmons, A., & Marsden, B. G. 1995, Icarus, 116, 
180 

Wurm, G., & Blum, J. 1998, Icarus, 132, 125 



This preprint was prepared with the AAS IATeX macros v4.0. 



34 



Table 1. Basic Model Parameters 



Parameter 


Symbol 


1 AU Models 


35 AU Models 


Width of annulus 


da 


0.17 AU 


6 AU 


Initial velocity 


V 


4.7 m s _1 


0.45-45 m s" 1 


Particle mass density 


Po 


3 g cm -3 


1.5 g cm -3 


Relative gas velocity 


T] 


60 m s _1 


30 m s" 1 


Time step 


St 


0.5 yr 


1-250 yr 


Number of mass bins 


N 


70-400 


64-256 


Mass spacing of bins 


6 


1.1-2.0 


1.25-2.0 


Minimum velocity for cratering 


Vf 


50 m s" 1 


1 cm s _1 


Impact strength 


So 


10 7 erg g _1 


2 x 10 6 erg g _1 


Crushing energy 


Qc 


10 8 erg g _1 


5 x 10 7 erg g _1 


Fraction of KE in ejecta 


Ike 


0.10 


0.05 


Coefficient of restitution 


Cl 


lO" 7 


lO" 2 


Coefficient of restitution 


C2 


lO" 7 


10" 3 
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Table 2. Model Results at 1 AU 



S = 1.12 5 = 1.25 5 = 1.4 ,5 = 2.0 



Time (yr) 


fmax 


Tfl ( V m a x ) 


N 




fFI> i^max ) 


N 


Tmax 


Tfl ( V m a x ) 


N 


Tmax 


( T m a x ) 


N 


1.0 x 10 3 


26.8 


2.4 xlO 20 


7 


25.7 


2.1 XlO 20 


5 


23.6 


1.6 


XlO 20 


33 


19.1 


8.7 


xlO 19 


646 


7.1 x 10 3 


326.0 


4.4 XlO 23 


2 


304.1 


3.5 XlO 23 


1 


284.1 


2.9 


xlO 23 


1 


153.0 


4.5 


xlO 22 


10 


1.0 x 10 4 


525.2 


1.8 XlO 24 


1 


509.3 


1.7 xlO 24 


2 


488.6 


1.5 


xlO 24 


1 


314.8 


3.9 


xlO 23 


8 


2.0 x 10 4 


952.2 


1.1 xlO 26 


3 


976.5 


1.2 xlO 26 


1 


873.2 


8.4 


xlO 24 


1 


801.5 


6.5 


xlO 24 


5 


6.3 X 10 4 


1450.4 


3.8 XlO 26 


1 


1470.2 


4.0 xlO 26 


2 


1446.0 


3.8 


xlO 26 


5 


1675.7 


6.1 


xlO 26 


1 


1.0 X 10 6 


1931.4 


9.1 XlO 26 


1 


1780.5 


7.1 xlO 26 


1 


1830.7 


7.7 


xlO 26 


3 


1915.0 


8.8 


xlO 26 


1 


1.2 X 10 6 


1962.7 


9.5 XlO 26 


1 


1801.4 


7.4 xlO 26 


2 


1888.2 


8.5 


XlO 26 


1 


1948.0 


9.3 


xlO 26 


1 



Notes to Table 2. 

These results are for the WS93 fragmentation algorithm. The columns list the mass and radius of the N largest bodies, m(r max ) in g and 
Tmax in km, as a function of time for four values of the mass spacing, 8. 
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Table 3. Model Results at 1 AU 



Time (yr) 




5 = 1.25 






S = 1.4 






5 = 2.0 




fmax 


WT> ij'max ) 


N 


Tmax 


WT> ij'max ) 


N 


fmax 


WT> ij'max ) 


N 


1.0 x 10 3 


25.5 


2.1 xlO 20 


7 


23.6 


1.6 xlO 20 


33 


19.1 


8.7 xlO 19 


646 


7.1 x 10 3 


283.1 


2.9 xlO 23 


5 


257.4 


2.1 xlO 23 


5 


159.3 


5.1 xlO 22 


44 


1.0 x 10 4 


523.3 


1.9 xlO 24 


1 


474.9 


1.3 xlO 24 


2 


364.4 


6.1 xlO 23 


1 


2.0 x 10 4 


948.1 


1.1 xlO 26 


1 


948.1 


1.1 xlO 26 


1 


812.5 


6.7 xlO 24 


2 


6.3 x 10 4 


2068.8 


1.1 xlO 26 


1 


1739.5 


6.6 xlO 26 


1 


1479.2 


4.1 xlO 26 


4 


1.0 x 10 6 


2316.6 


1.6 xlO 20 


1 


2129.6 


1.0 xlO 26 


1 


1724.8 


6.5 xlO 26 


7 


1.2 x 10 6 


2407.3 


1.8 xlO 20 


1 


2212.6 


1.2 xlO 26 


1 


1793.9 


7.3 xlO 26 


7 



Notes to Table 3. 

These results are for the Davis et al. fragmentation algorithm. The columns list the mass and radius of the N largest bodies, 
m(rmav) in g and r ma x in km, as a function of time for three values of the mass spacing, S. 
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Table 4. Fragmentation Results at 100 Myr for e = 10 3 (No Velocity Evolution) 



M (M E ) 


CR 


r (m) 


r 9 5% (km) 


r 5 (km) 


r max (km) 


r (j 'max ) 


T P (Myr) 


1 


IO" 5 


80 


0.3 


0.6 


0.6 


2660 


2764 


3 


IO" 5 


80 


1.0 


1.6 


1.8 


46 


919 


10 


10~ B 


80 


3.7 


5.8 


7.0 


2034 


276 


30 


io~ B 


80 


657.2 


337.2 


1006.0 


3 


92 


100 


io~ B 


80 


587.8 


451.7 


1022.5 


1 


27.6 


1 


io- 2 


80 


0.3 


0.6 


0.6 


1223 


2803 


3 


10- 2 


80 


1.0 


1.6 


1.9 


2 


937 


10 


IO" 2 


80 


3.5 


5.8 


6.6 


61 


280 


30 


IO" 2 


80 


659.3 


333.9 


1080.1 


1 


93.7 


100 


IO" 2 


80 


659.4 


477.6 


1133.4 


2 


27.6 



Notes to Table 4. 

The columns list rg 5 % (95% of the mass is contained in objects with radius less than r a5 %), r$ (the radius where the 
cumulative number of objects is 10 B ), r max (the radius of the largest object), N(r max ) (the number of objects with 
radius r max , and Tp (the timescale to produce 1000 km objects) as functions of the initial mass in the annulus (Mo), the 
coefficient of restitution (cr), and the initial radius of the largest object (ro). 



38 



Table 5. Fragmentation Results at 100 Myr for e = 10 3 (Limited Velocity Evolution) 



M (M E ) 


CR 


r (m) 


r 9 5% (km) 


r 5 (km) 


r max (km) 


r (j 'max ) 


T P (Myr) 


1 


lO" 2 


80 


0.3 


0.6 


0.6 


1084 


2150 


3 


10- 2 


80 


1.0 


1.6 


1.9 


3 


720 


10 


lO" 2 


80 


3.9 


6.4 


7.5 


31 


216 


30 


lO" 2 


80 


656.1 


333.2 


1089.9 


1 


72.1 


100 


lO" 2 


80 


735.0 


537.4 


1240.4 


4 


21.6 



Notes to Table 5. 

The columns list r 95 % (95% of the mass is contained in objects with radius less than r g5 %), r$ (the radius where the 
cumulative number of objects is 10 5 ), r max (the radius of the largest object), N(r max ) (the number of objects with 
radius r max , and Tp (the timescale to produce 1000 km objects) as functions of the initial mass in the annulus (Mo), the 
coefficient of restitution (cr), and the initial radius of the largest object (ro). 
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Table 6. Fragmentation Results for c R = 10 2 (Complete Velocity Evolution) 



Mo (M E ) 


eo 


r (m) 


r 95% (km) 


r B (km) 


r max (km) 


N(r max ) 


t p (Myr) 


1 


io~ 4 


80 


1.6 


22.1 


1000.1 


1 


447.8 


10 


10" 4 


80 


0.8 


49.8 


1000.2 


4 


20.1 


1 


IO" 3 


80 


3.5 


26.0 


997.2 




892.7 


3 


io- 3 


80 


2.5 


32.7 


1000.2 




183.6 


10 


io- 3 


80 


1.8 


52.6 


1036.8 




36.5 


30 


IO" 3 


80 


1.8 


73.0 


1002.6 




10.4 


100 


IO" 3 


80 


1.8 


115.1 


1143.1 




2.6 


10 


10- 2 


80 


17.7 


50.4 


1000.2 


1 


428.3 


100 


io- 2 


80 


11.3 


89.0 


1004.3 


5 


25.4 



Notes to Table 6. 

The columns list r 95 % (95% of the mass is contained in objects with radius less than r g5 %), r$ (the radius where the 
cumulative number of objects is 10 5 ), r max (the radius of the largest object), N(r max ) (the number of objects with radius 
fmax), and Tp (the timescale to produce 1000 km objects) as functions of the initial mass in the annulus (Mo), the 
coefficient of restitution (cr), and the initial radius of the largest object (ro). 
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Table 7. Fragmentation Results as a Function of Strength and Mass Spacing 



M [Me) 


eo 





So (erg g ) 


T f (Myr) 


r 5 (km) 




10 


10" 4 


2.0 


IO 1 


50 


55.4 


468.2 


10 


10" 4 


2.0 


IO 1 


100 


55.4 


473.4 


10 


10~ 4 


2.0 


10 2 


50 


54.7 


670.6 


10 


10~ 4 


2.0 


10 2 


100 


54.7 


776.1 


10 


IO -4 


2.0 


10 3 


50 


55.3 


1200.8 


10 


IO" 4 


2.0 


10 3 


100 


55.3 


1210.6 


10 


IO" 4 


2.0 


10 4 


50 


55.8 


1584.8 


10 


10~ 4 


2.0 


10 4 


100 


55.8 


1605.0 


10 


10~ 4 


2.0 


10 5 


50 


62.8 


2002.6 


10 


IO -4 


2.0 


10 5 


100 


63.2 


2331.3 


10 


IO -4 


1.4 


10 4 


50 


50.1 


1506.4 


10 


10~ 4 


1.4 


10 4 


100 


50.1 


1706.2 


10 


10~ 4 


1.4 


10 5 


50 


50.1 


1525.4 


10 


10~ 4 


1.4 


10 5 


100 


50.2 


2254.8 


10 


IO" 4 


1.4 


2 x 10 6 


50 


50.2 


1552.6 


10 


IO -4 


1.4 


2 x 10 6 


100 


50.5 


2382.8 


10 


io- 3 


2.0 


10 4 


50 


47.0 


1621.0 


10 


io- 3 


2.0 


10 4 


100 


47.1 


2038.2 


10 


IO" 3 


2.0 


2 x 10 6 


50 


52.8 


1452.9 


10 


IO" 3 


2.0 


2 x 10 6 


100 


53.2 


2019.2 


10 


io- 3 


1.4 


10 4 


50 


55.5 


1381.4 


10 


io- 3 


1.4 


10 4 


100 


55.6 


1834.9 


10 


io- 3 


1.4 


2 x 10 6 


50 


52.8 


1452.9 


10 


IO" 3 


1.4 


2 x 10 6 


100 


53.2 


2019.2 


10 


IO" 3 


1.25 


2 x 10 6 


50 


45.6 


2309.1 


10 


io- 3 


1.25 


2 x 10 6 


100 


45.6 


3333.0 



Notes to Table 7. 

The columns list (the radius where the cumulative number of objects is 10 s ) and r max (the radius of the largest object), 
N(r max ) (the number of objects with radius r max at two times, r/, as functions of the initial mass in the annulus (Mo), 
the initial eccentricity (eo) the mass spacing (<5) and the intrinsic tensile strength (So). 
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10 15 20 25 

log Mass (g) 



Fig. 1. — Cumulative mass distribution at selected times for a calculation at 1 AU with M 
= 0.667 Me and 5 = 1.25. A group of runaway bodies with mi ^ 10 24 g forms at ~ 10 4 yr. 
These batches contain 13% of the initial mass at 6.3 x 10 4 yr and 17% at 1.15 x 10 5 yr. 
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Fig. 2. — Velocity evolution for the 1 AU model in Fig. 1. (a) Horizontal velocity 
distribution. Viscous stirring increases all velocities with time; dynamical friction brakes 
the runaway bodies and increases velocities of the lowest mass bodies, (b) Ratio of vertical 
to horizontal velocity. The ratio remains close to the equilibrium value of V/H = 0.53 for 
all but the most massive objects, where V/H ~ 0.3-0.5. 
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Fig. 3. — Cumulative size distributions for a Kuiper Belt model with M = 10 M E and C\ 
= c 2 = 1CT 5 . The eccentricity is constant in time at e = 1CT 3 . Collisional growth is slow for 
200 Myr until the largest bodies have r max = 50 km. Runaway growth begins when r max ^ 
100 km. These bodies then grow to sizes of 10 3 km to 10 4 km in 20-30 Myr. The size 
distribution consists of a fragmentation tail with N c oc r~ 2 25 for r\ ^ 0.1 km, a transition 
region with a flatter size distribution, and an accretion component with Nc oc r^ 2 ' 6 . 
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Fig. 4. — Cumulative size distributions for a Kuiper belt model with M = 10 M E , c\ = 
1CT 2 , and c 2 = 1CT 3 . The eccentricity is constant in time at e = 1CT 3 . The evolution is 
identical to that in Figure 3, except for the appearance of a 'kink' in the size distribution 
due to rebounds at ~ 5 m. As the evolution proceeds, debris from collisions with ^> 
5 m suppresses this kink and the size distribution more closely follows the N c oc r~ 2 ' 5 derived 
in calculations without rebounds. 
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Fig. 5. — Cumulative size distributions for a Kuiper belt model with M = 10 M E , c\ = 1CT 2 , 
C2 = 10~ 3 , and limited velocity evolution. Eccentricity changes are due to fragmentation, 
which places a fixed fraction f KE of the kinetic energy of each collision into the debris. The 
size distribution develops a kink at r\ m 5 m due to rebounds at 10-50 Myr when r max ps 
1-3 km. As the maximum radius increases, the size distribution approaches a smooth power 
law distribution, N c oc r^ 2 ' 5 from r\ = 1 m up to r\ = 300-500 km. 
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Fig. 6. — Horizontal velocity distributions as a function of time for the limited velocity 
evolution model in Figure 5. The velocity dispersion of all bodies begins at hi — 4 m s _1 . As 
the maximum radius increases from r max = 80 m to r max = 1 km, collision debris increases 
the velocity of low mass bodies. This increase is most pronounced for bodies with r\ ~ 10 
m at early times, because rebounds reduce the velocities of the lowest mass objects (dotted 
line). The high velocity dispersion propagates to bodies with = 1-100 m as rebounds 
become less frequent and as mergers of the largest bodies produce more debris (dot-dashed 
and dashed lines). Once r max m 1000 km, the velocity distribution has settled into a high 
velocity component at ~ 1-100 m, a transition component at ~ 0.1-3 km, and a group 
of large objects, r\ ^ 3-5 km with hi ~ 3 m s -1 . 
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12 3 

log Time (Myr) 

Fig. 7. — Maximum radius as a function of initial mass for models with no velocity evolution 
(thin solid lines) and limited velocity evolution (thick solid lines). The initial eccentricity 
is e = 1CT 3 . The time to produce 1000 km objects scales inversely with initial mass, 
t p ~ r (M /10M E ) _1 , with r xs 276-280 Myr for models with no velocity evolution and 
t 216 Myr for models with limited velocity evolution. Fragmentation speeds up the 
growth of models with limited velocity evolution by redistributing kinetic energy from large 
objects to small objects. 
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Fig. 8. — Size and velocity distributions for a model with M = 10 Me, e = 10 -3 , 
= 2 x 10 6 erg g _1 , and velocity evolution: (a) cumulative size distribution, and (b) 
horizontal velocity as a function of time. Collisional growth is slow until the largest bodies 
have r max = 1-2 km at 9-10 Myr. Collisional damping reduces the velocities of all bodies to 
~ 1-2 m s _1 on this timescale; dynamical friction reduces the velocities of larger bodies to 
~ 10~ 2 m s _1 . Runaway growth then produces objects with radii of 100 km in another 2-3 
Myr. Viscous stirring increases particle velocities as objects grow to sizes of 100-300 km, 
and runaway growth ends. A prolonged linear growth phase leads to the production of 1000 
km objects; the horizontal velocities are then ~ 30-40 m s -1 for the smallest objects and ~ 
1 m s _1 for the largest objects. 
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Fig. 9. — Evolution of the maximum radius, r max , with time as a function of initial mass, 
M , for e = 1CT 3 . The timescale to reach runaway growth decreases with increasing M . 
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Fig. 10. — Evolution of r 5 , the radius where the cumulative number of objects is 10 5 , with 
time as a function of initial mass, M , for eo = 1CT 3 . The constraint on r 5 set by current 
observations is indicated by the horizontal dashed line. 
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Fig. 11. — Cumulative size distribution when r max = 1000 km for three values of the initial 
eccentricity, eo- Each size distribution follows two power laws, Nc oc r -2 ' 5 at small radii and 
Nc oc r~ 3 at large radii. The transition between the two power laws moves to larger radii 
as e increases. 
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Fig. 12. — Maximum radius as a function of strength for models with Mq = 10 Me and 
various eo and 5 as listed in the legend. 
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Fig. 13. — Cumulative size distributions for a model with Mq = 10 Me, eo = 10~ 4 , and So = 
100 erg g -1 . Runaway growth produces objects with r max m 100 km as in the standard model 
with S = 3 x 10 6 erg g -1 . Catastrophic disruption begins when r max ps 400 km. Continued 
mass loss due to catastrophic disruption reduces the population of low mass objects with 
Ti 0.1 km, which are completely removed from the distribution for r ^ 100 Myr. 
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